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1 Abstract 

The accuracy is studied of various numeri- 
cal flux functions for the inviscid fluxes when 
used for Navier-Stokes computations. The 
flux functions are benchmarked for solutions 
of the viscous, hypersonic flow past a 10° 
cone at zero angle of attack using first-order, 
upwind spatial differencing. The Ilarten- 
Lax/Roe flux is found to give a good bound- 
ary layer representation, although its robust- 
ness is an issue. Some hybrid flux formulas, 
where the concepts of flux-vector and flux- 
difference splitting are combined, are shown 
to give unsatisfactory pressure distributions; 
there is still room for improvement. Investi- 
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gations of low diffusion, pure flux- vector split- 
tings indicate that a pure flux-vector split- 
ting can be developed that eliminates spuri- 
ous diffusion across the boundary layer. The 
resulting first-order scheme is marginally sta- 
ble and not monotone. 


2 Introduction 

This paper addresses the problem of design- 
ing accurate numerical flux functions ap- 
proximating the inviscid fluxes in the Euler 
and Navier-Stokes equations, and may be re- 
garded as a sequel to [1]. In the latter paper it 
was demonstrated, among other things, that 
numerical flux functions that do not recog- 
nize contact and shear discontinuities grossly 
exaggerate the diffusion of entropy and shear 
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across the boundary layer in a Navier-Stokes 
calculation. Specifically, flux formulas of the 
flux-vector splitting (FVS) type, such as Van 
Leer’s [2], and those based on a tuned scalar 
viscosity coefficient, such as Jameson’s [3], 
were found to be inferior to those based on 
flux-difference splitting (FDS), such as Roe’s 

[4] - 

Since the appearance of [1], a clear shift in 
the use of the various flux formulas has been 
observed. FVS is gradually being phased out 
as a component of Navier-Stokes codes in fa- 
vor of FDS, while the scalar viscosity coef- 
ficient in central-difference schemes on occa- 
sion, [5], has been replaced by a viscosity ma- 
trix , yielding the same low diffusion as FDS. 

Another trend, inspired by developments 
in hypersonic flight, has been the extension 
of known flux formulas for ideal gases to real 
and reacting gases. Examples can be found 
in [6, 7, 8]. 

At the same time, a growing effort is being 
spent on the development of genuinely multi- 
dimensional schemes, in which the influence 
of the grid coordinate directions is reduced 
as much as possible. This brings along the 
formulation of multi-dimensional flux func- 
tions, such as the one due to Rumscy et al. 
[9, 10], which is based on an approximate Rie- 
mann solver including waves traveling in two 
orthogonal directions of physical importance. 
Another example is the one by Goorjian and 
Obayashi [11], based on waves traveling in 
and normal to the flow direction. 

Independent of the above developments, 
several attempts have still been made to sal- 
vage the concept of FVS, mainty because of 
three reasons: 

1. the formulas are relatively simple; 

2. the split fluxes are easy to linearize, for 
the benefit of implicit marching schemes, 

3. the extension to real gases is relatively 
straightforward 


Worth mentioning in this field is the work 
of Hanel [12, 13], who has suggested a num- 
ber of modifications to Van Leer’s [2] FVS, 
including a mix with FDS formulas. This 
hybrid approach was carried further by Van 
Leer [14] and appears to be completed by 
Liou and Steffen [15], whose AUSM (Advec- 
tive Upwind Splitting Method) flux appears 
to rival the accuracy and robustness of Roe’s 
at significantly reduced computational com- 
plexity. 

In the present paper we review some of 
these new developments in FVS, and subse- 
quently pose the following question: 

• Is it possible at all to construct a pure 
FVS that does not diffuse a grid-aligned 
boundary layer , and makes a stable 
combination with some form of time- 
marching? 

The answer turns out to be “yes,” but loss 
of monotonicity of the numerical solution is 
unavoidable. 

Thus, improvements to date of the origi- 
nal Van Leer FVS can be grouped into two 
categories: 

• ’’Hybrid”, or mixed FVS/FDS modifica- 
tions 

• ’’Pure” FVS modifications 

To compare the accuracy of the various 
flux functions when used for a viscous com- 
putation, the viscous, hypersonic flow over 
a cone at zero angle of attack is used to 
benchmark the flux formulas. In keeping 
with the style of the earlier paper, [1], we 
present solutions for the Mach 7.95 viscous 
flow over a 10° cone at a Reynolds number 
of Rei = 0.42 x 10 6 and freestream total 
temperature of T 0?oo = 775.56AL Adiabatic 
wall boundary conditions are applied at the 
cone surface, resulting in a wall temperature 
of T w — 11.73 Tqo and a boundary layer thick- 
ness of approximately 0.5°. In contrast to the 


2 


approach in [1], which used a two-dimensional 
code for the flux comparisons, we solve the 
one-dimensional conical Navier-Stokes equa- 
tions. In all of the following computations, 
the residual is formed using first-order up- 
wind differencing with the various flux func- 
tions. Unless otherwise noted, Yoon’s [16] 
LU-SGS scheme for approximate Newton it- 
eration is used to obtain the solutions. All of 
the computations reported here were made on 
a uniform mesh composed of 50 cells spanning 
5° from the cone surface. For the benchmark 
cases, this results in approximately 8 points 
in the boundary layer. The Harten-Lax/Roe 
flux that was presented, but not tested, in [1] 
is compared to Roe’s FDS scheme for the con- 
ical flow in Section 3. Section 4 outlines re- 
cent hybrid modifications to the original FVS 
scheme and compares these to the original 
FVS and Roe’s FDS scheme for the conical 
flowfield. Section 5 addresses the prospect 
of pure FVS splitting and shows that, al- 
though the adverse dissipative properties of 
the original scheme can be negated, the re- 
sulting flux formulas are non-monotone. This 
non-monotonicity is shown to be unavoidable, 
and is evident from an examination of the 
eigenvalues of the split flux Jacobians as well 
as the eigenvalues of a linearized representa- 
tion of the residual. 

3 Harten-Lax/Roe Flux 

The Harten-Lax/Roe flux ([17] and [18]) is an 
FDS formula, incorporating a “smart” scalar 
dissipation coefficient: the scalar is a square- 
amplitude weighted average of the character- 
istic speeds. This formula was discussed in 
[1], but not tested. The Harten-Lax/Roe flux 
does not yield linear stability, as the scalar 
dissipation coefficient may be too small for 
stabilizing the weaker waves with the larger 
characteristic speeds. The instability will 
show up first in these waves, increasing their 


amplitude, which then feeds back into the dis- 
sipation coefficient, increasing its value; thus, 
stability is restored. 

Numerical solutions for the temperature 
and the pressure are shown in Figures 1 and 
2; for comparison the results of our bench- 
mark flux, i.e. Roe’s [4], are also given. 
As can be seen from the figures, agreement 
with the benchmark results is excellent, al- 
though very minor pressure oscillations are 
observed in the boundary layer. The pres- 
ence of these oscillations might indicate a ten- 
dency towards a non-linear instability that 
could show up in a more demanding calcula- 
tion. Convergence was readily achieved with 
the LU-SGS scheme using an infinite Courant 
number. Convergence was also achieved us- 
ing a single stage explicit (forward Euler) 
scheme, but a very small Courant number 
was needed. This convergence behavior was 
found on a range of 32-bit machines, while 
on two 64-bit machines (CRAY Y-MP and 
X-MP), this flux combined with either the 
LU-SGS or forward Euler scheme yielded only 
a three order of magnitude drop in the L 2 
norm of the residual. The Harten-Lax/Roe 
flux nominally uses all of the wave strengths 
and speeds evaluated at Roe’s averaged state, 
but blends these into a single, weighted wave 
speed for the construction of the numerical 
flux, ft is only a small step computing Roe’s 
flux once the wave strengths and speeds are 
found, and one might as well go all the way, 
and compute the flux using Roe’s formula. As 
pointed out in the previous paper [1], though, 
the weighted speed can be obtained more sim- 
ply as a ratio of scalar products: 

AWo AF 

V = (1) 

AWoAU 

where U is the vector of conserved state 
quantities, F is the flux vector, and W is an 
alternative state vector, viz. the gradient of 
an entropy function associated with the Euler 
equations; for instance W = ■ In spite 
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Figure 1: Conical Navier-Stokes solution us- Figure 2: Same solution as in Figure 1; pres- 
ing Harten/Lax/Roe’s scheme; temperature sure distributions, 
distributions. 


of this simplification, this ingenious scheme 
might not obtain the notoriety that Roe’s 
scheme has obtained because of the stability 
hazard. 

4 Hybrid Modifications 
of Van Leer’s FVS 

We recall the definition of FVS: 

F(U) =F+(U) + F-(U); (2) 

here U is the vector of conserved state quan- 
tities, F(U) is the vector of inviscid fluxes 
in one coordinate direction, and F + (U) and 
F”(U) are called the forward and backward 
fluxes, respectively. 

The FVS used most frequently in practice 
is Van Leer’s [2], owing to the following design 
features: 

1. The split fluxes are continuously differ- 
entiable, which preserves the numerical 
accuracy near sonic points, and allows 
smooth linearization; 


2. For subsonic flow, the 

Jacobians dF + /dU and cIF~ /dXJ have a 
zero eigenvalue, which accounts for crisp 
numerical profiles of steady shocks. 

If the second constraint is relaxed, a one- 
parameter family of continuously differen- 
tiable split fluxes can be generated [6]; these 
are the simplest possible in the sense that 
they are at most quartic in the Mach number, 
just as the Van Leer fluxes. The differences 
among members of this family arise only in 
the energy-flux splitting. 

Included in this family is the energy-flux 
splitting originally proposed by Ilanel [12]: 


where H is the specific total enthalpy. Ad- 
vantages of this formula are: 

1. it is as simple as can be; 

2. it admits steady Euler solutions wjth 
constant total enthalpy throughout the 
flow. 

Hand claims that this flux splitting, when 
used in Navier-Stokes calculations, gives 
more accurate total-enthalpy values in the 
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boundary layer. This may have been ob- 
served for the lower flow speeds; in the hy- 
personic flow regime the improvement is in- 
significant. 

The problem of numerical diffusion across 
a boundary layer, even when it is aligned 
with the grid, can easily be understood by 
rewriting the formulas for the flux of a general 
scalar function, <f>. For the inviscid fluxes, <f) 
would be replaced by u and H for the convec- 
tive transverse momentum and energy fluxes, 
respectively. Letting D denote the bottom, 
and T denote the top cell: 

{p v )<f> = + F- T <f>T 

= (F*, e + F-, T )i^L 

-(F+b - F-, t )^^L ( 4 ) 

The first term in these expressions represents 
a central-differencing flux; numerical diffu- 
sion is introduced by the second term. When 
the net mass flux 

FT = et.B + F-J (5) 

vanishes, this is because of a cancellation , not 
because F^ B and F^ T vanish individually. 
As a consequence, the diffusive terms do not 
vanish with the mass flux. 

Hanel [13] has suggested to replace the for- 
mula for the transverse-momentum flux by 
one borrowed from flux-difference splitting: 


(pu)u — F m u U p W i n( j, 

> 

( 6 ) 

^upwind — ^ B 

if F net 

> 0 , 

( 7 ) 

^upwind ~ = ^ T 

if F net 

< 0 . 

( 8 ) 


This mixture of flux- vector splitting and flux- 
difference splitting prevents the numerical 
broadening of the boundary layer, but can 


not improve the accuracy of the wall temper- 
ature. It further introduces pressure irregu- 
larities across the boundary layer. 

The similarity of the transverse momen- 
tum and energy fluxes, though, suggests that 
a further improvement can be obtained by 
introducing a similar formula for the energy 
flux, i.e.: 

(pv)H = p, 4 „ d , (9) 

with 

^upwind = Ih if F m l — 0, (10) 

Upwind = H t ifF" et <0. (11) 

As reported in [14], this removes the error 
in the wall temperature; unfortunately, the 
pressure irregularity remains. 

The next pair of figures shows the tempera- 
ture and pressure distributions obtained with 
the three variations of Van Leer’s FVS repre- 
sented by Equations (3), (6) and (9); for com- 
parison, results for the original Van Leer FVS 
and for Roe’s FDS are included. Use of the 
first modification of the energy flux, Hanel ’s 
energy flux, (equation (3), ’’Hanel 87” in the 
Figure), hardly causes any change with re- 
spect to using the Van Leer flux; the second 
modification, Hanel’s FDS-like transverse- 
momentum flux, (equation (6), ’TTanel 89” in 
the Figure), does achieve the proper shrink- 
ing of the boundary layer, which also corrects 
the pressure profile (except for a fluctuation 
near the wall), but the wall temperature is 
still wrong. The third modification, Van Leer 
and Hanel ’s FDS treatment of both the trans- 
verse momentum and energy fluxes, (equa- 
tions (6) and (9), ”VL/H 90” in the Figure), 
finally corrects the wall temperature, but the 
pressure fluctuation remains. 

The final stage in this sequence of hybrid 
flux formulas is the technique named AUSM, 
developed by Liou and Steffen [15], where all 
advective terms in the fluxes are treated using 
the advective Mach number splitting from the 


5 




original Van Leer mass splitting. The remain- 
ing pressure term is split as usual. Prelimi- 
nary numerical testing on the cone-flow prob- 
lem indicates that the results are comparable 
Hanei F 87 S to those of Roe’s flux; in particular, they show 
vl^h 90 a smooth pressure distribution. Some other 

Roe r DS 

results are included in the forum-paper by 
Liou and StefFcn [19], included in the present 
conference proceedings. Since this scheme 
was not available in the literature at the time 
that this paper was written, it will not be 
included in the comparisons made here. 


5 Pure FVS Modifica- 
tions 


Figure 3: Conical Navier-Stokes solutions us- 
ing hybrid modifications to the original Van Unaffected by the development of hybrid 
Leer FVS; temperature distributions. FVS/FDS formulas, an intriguing question 

still remains: 


• Is it possible at all to construct a 
pure FVS that does not diffuse a grid- 
aligned boundary layer and makes a sta- 
ble combination with some form of time- 
marching? 



Orig. FVS 
Hand 87 
Hanel 89 
VL/H 90 
Roe FDS 


This question can be rephrased as: 

• Is it possible to split the Elder fluxes such 
that both F+ and F t ~ vanish with the flow 
speed } while numerical stability is main- 
tained? 

If indeed this were possible, the form of the 
flux splitting for small v would follow imme- 
diately from symmetry considerations: 


Figure 4: Same solutions as in Figure 3; pres- 
sure distributions. 


II 

l -pv + 0(d), 

(12) 

F ± — 

x mom|| 

tp ± spav + 0(d), 

(13) 

— 

moml 

l -puv + 0(d), 

(14) 

F ± = 

A energy 

\pHv + 0(d), 

(15) 
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where a is the sound speed and s is a free 
parameter representing the derivative of the 
split pressures for v = 0. Most noticeable 
is that for v = 0 this splitting leads to cen- 
tral differencing, which will be unstable if 
forward time-differencing is used. If s also 
vanishes, central differencing will spread to a 
small neighborhood of v = 0. Positive values 
of s would seem to introduce some dissipa- 
tion, since this will make the split pressures 
upwind-biased; this turns out to be not nec- 
essarily true. The fluxes (12-15), valid for 
v « 0, must smoothly join the branches for 
larger values of \v\; for the latter we may use 
the standard Van Leer fluxes. 

Below we shall study a four-parameter fam- 
ily of splittings of the one-dimensional Euler 
fluxes, hereafter referred to as Low Diffusion 
FVS (LDFVS), defined by the formulas, valid 
for |M| < 1 


2^ 
3 H- 

= ±/>a(t(M± l) 2 



1 

t\3 

1 

(16) 

F L m 

= F^v + p*, 

(17) 

p ± 




T >(M 2 - If) , 

(18) 

F ± 

energy 




^fupa 3 M 2 (M‘ 2 - l) 2 

(19) 


The three chief parameters are //, v and u>. 
If/i = i/ = u; = 0, the formulas return the 
standard Van Leer fluxes, with Hanel’s en- 
ergy flux, regardless of the value of P. The 
higher the value of the exponent P, the nar- 
rower the interval around M = 0 where the 
fluxes deviate significantly from the Van Leer- 
Hanel fluxes. For conciseness, we will only 
show results for P = 2 since the results for 
higher P are similar in nature. Note that the 
momentum-splitting is conceived as a split- 



Figure 5: Normalized forward mass flux, 

F+/(/kz), for three values of //, with P — 2. 


ting of the convective flux, plus a pressure 
splitting. The extra terms in the energy split- 
ting are considered less important than those 
in the mass and momentum fluxes, since the 
splitting should also be valid for isothermal 
flow (7 = 1); in this case the energy equation 
would drop out completely. Similarily, inclu- 
sion of a transverse-momentum flux is not a 
top priority, since the splittings should be sta- 
ble in the first place for a one-dimensional 
flow. 

Figures 5, 6 and 7 show the mass, pres- 
sure and energy splittings of this family. In 
Figure 5, three values of the mass splitting 
parameter, /i, are shown for the forward flux, 
where p = 0 results in the original Van Leer 
splitting and // = 1 causes the split mass flux 
to be identically zero at zero Mach number. 
Note that for the larger values of // the split 
flux becomes negative for a range of negative 
Mach numbers, which is somewhat odd for a 
’’forward” flux, and could be destablizing. 

Figure 6 shows the split pressure from the 
forward momentum flux. For v ~ 7 the 
slope of p + vanishes at M — 0, correspond- 
ing to s = 0 in Equations (12-15); as men- 
tioned above this leads to central pressure- 
differencing in the neighborhood of M = 0. 
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1.00 



v = 1 .5 


- 1.00 


Figure 6: Normalized “forward” pressure 

for three values of p, with P = 2. 


This value unexpectedly performed well re- 
garding stability. For v > | the split pres- 
sures lose monotonicity, which has little phys- 
ical appeal; that this is an improper choice is 
~ *v = 3/4 a l so borne by a stability analysis. 

Figure 7 shows the resulting energy split- 
tings for the case of vanishing split mass 
fluxes at M = 0 (he. // = 1), for selected 
values of the energy-splitting parameter, to. 

Liou and Steffen in [20] have independently 
derived a pure FVS formula along the same 
lines as presented above. Their flux has been 
coined HOPE (for ITigh-Order Polynomial 
Expansion) and includes splittings that are 
similar to those presented here. The mass, 
pressure and energy splittings of the HOPE 
family are 






±^{(M±lf 

+m,(M 2 - l) 2 } 
F^v + p*, 

2 '1 


( 20 ) 

( 21 ) 


pa 

7 


-(M ±\) 2 {2t M) 

4 




if}, 



F ± 

energy 


= F* H 

mass 


( 22 ) 

(23) 


Figure 7: Normalized forward energy flux, 
/(pa 3 ), for three values of u 


energy 


where 

mj(M) = (M 2 - 1)/(M 2 + l) 5 (24) 

A comparison with Eqs. (16-19) shows that 
the HOPE splitting is very similar to the LD- 
FVS splitting with p = 1, v = | and u> = 0. 
The function ?ni(M) is a blending function 
that ensures that the split mass flux is zero 
at M = 0 and that the split fluxes smoothly 
join the unsplit fluxes near \M\ — 1. As sug- 
gested in [20] the parameter S in (24) is taken 
to be S — 4. In [20], various pressure split- 
tings were tested, where the one shown here 
was found to be the most robust. 

Figures 8 and 9 compare mass and pres- 
sure splittings of the HOPE and LDFVS fam- 
ilies. The effect of the blending function 
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Figure 8: Normalized forward mass flux for 
the HOPE and LDFVS (/< = 1) splittings. 



Figure 9: Normalized ” forward” pressure flux 
for the HOPE and LDFVS (v = f splittings. 


in HOPE is apparent in the Figures. The 
pressure splittings are much the same, with 
the HOPE curve showing a narrower plateau 
around M — 0 and closer agreement with the 
Van Leer curve near |A/| = 1. The HOPE for- 
ward mass flux follows the unsplit flux closely 
for M > 0 and oscillates about F+ = 0 for 
M < 0, while the LDFVS scheme is overall 
more smooth. Both pressure and mass split- 
tings of the LDFVS family can be brought 
close to those of the HOPE family ( S = 4) 
by using a higher exponent P, viz. P = 6. 
As will be shown in a following section, the 
smoothness (or lack of smoothness) of the 
split fluxes has a direct influence upon the 
smoothness of the split-flux eigenvalues. 

The similarity of the above flux splittings 
yields comparable stability bounds and, as 
will be shown in the next section, both 
fluxes can lead to non-monotone solutions, 
even with first-order upwind differencing. 
This non-monotonicUy is not readily appar- 
ent from the conical flow results shown be- 
low, owing in part to the geometric and vis- 
cous source terms appearing in the equa- 
tions, but it appears to be unavoidable in 
multi-dimensional calculations. Regardless of 
this unsettling finding, the conical-flow re- 
sults show that the excess diffusion is elim- 
inated with these flux splittings. 

Results obtained with first-order upwind 
differencing incorporating the LDFVS fluxes 
(16), (17), (IS) and (19) are presented in Fig- 
ures 10 and 11; results for Van Leer and 
Roe fluxes again arc supplied for comparison. 
The splitting parameters used were: ft — 1 , 
v — to = 0 and P = 2. It will be shown 
in the following section that this combina- 
tion lies outside the stability range predicted 
by a stability analysis for inviscid flow. A 
converged solution was nevertheless obtained 
with the LU-SGS implicit marching scheme, 
while convergence could not be achieved with 
a single stage (forward Euler) explicit scheme. 
In fact, the results are the best obtained so 
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Figure 10: Conical Navier-Stokes solutions 

using the LDFVS formula for // = 1, v = |, Figure 12: Error in computed wall tempera- 
te = 0, and P = 2; temperature distributions, ture using various levels of fi for the conical 

Navier-Stokes solution. 



Figure 11: Same solutions as in Figure 10; 
pressure distributions. 

far with a first-order scheme. Presumably 
this is so because the scheme reverts to cen- 
tral differencing wherever v is small, i.e. in 
the boundary layer. The wall temperature 
is essentially the same as for the Roe flux; 
the pressure curve actually is better than for 
Roe’s flux, as it lacks the pressure dip appar- 
ent in the Roe results just above the bound- 
ary layer. 

To illustrate the effect of fi upon the dif- 
fusive error in the solution, we then ran the 
same test problem for 11 values of ft com- 
bined with the optimum v - value in the pres- 


sure splitting as found from a stability anal- 
ysis (Section 6), and lo = 0. All of the 
computations were made with the LU-SGS 
scheme, and convergence to machine zero was 
obtained in 750 or fewer iterations. In Figure 
12 the relative error in the wall temperature 
is plotted against //, where the baseline wall 
temperature is found from 

Tw.'x = Tc „ (l + sfPr^Ml) (25) 

The effect of the artificial diffusion is clear: 
for the uniform grid used, with about 8 cells 
in the boundary layer, the diffusion coefficient 
must be cut at least by an order of magnitude, 
i.e. 0.9 < // < 1, in order to admit a wall 
temperature of acceptable accuracy. 

Results obtained with the HOPE splitting, 
Equations (20), (21), (22) and (23), for the 
conical flow are presented and compared to 
other results in Figures 13 and 14. Since the 
mass, pressure and energy splittings are sim- 
ilar to those in the new FVS scheme, it is not 
surprising that the computed solutions are 
also very similar. Again, the LU-SGS scheme 
was used to converge the solution, and con- 
vergence to machine zero was obtained in ap- 
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Figure 13: Conical Navier-Stokes solutions Figure 14: Same solution as Figure 13; Pres- 
using the HOPE and LDFVS splittings: tem- sure distributions, 
perature distributions. 


proximately 1000 iterations. Here too conver- 
gence could not be obtained using an explicit 
time-marching scheme. This stability behav- 
ior is also indicated by the stability analysis 
below. 

It would appear from the above results for 
the one-dimensional conical flow that the new 
flux functions have potential for use in more 
realistic, multi-dimensional problems. This, 
however, turns out not to be true. To test the 
pure FVS schemes for a two-dimensional flow, 
the inviscid transonic flow about a NACA 
0012 at 1.5° angle of attack and free-stream 
Mach number of M = 0.85 was computed on 
a relatively coarse (64 by 32 cell) grid. It 
was impossible to obtain any solution with 
either the LDFVS splitting (// = 1) or the 
HOPE splitting, regardless of Courant num- 
ber or level of implicit dissipation (controlled 
by the splitting parameter in the LU-SGS 
scheme). For // = 0.9 a converged solution 
was obtained, although the resulting flowfield 
was highly non-monotone. To investigate this 
further, both fluxes were used to compute the 
M = 2 inviscid flow over a 10° compression 
ramp on a 50 by 50 cell grid. Both computa- 


tions converged to machine zero, but yielded 
highly non-monotone results downstream of 
the ramp shock. This finding is discourag- 
ing, but not surprising, as the stability and 
monotonicity analysis in the following section 
will show. Note that this non-monotonicity 
occurred with first-order spatial differencing 
and can not be cured by the kind of limiters 
used in higher order schemes. 

6 Stability and Mono- 
tonicity Analysis 

The usual stability argument regarding FVS 
is that B+ = dF+/d\J and B~ = d¥~/dU 
must have non-negative and non-positive 
eigenvalues, respectively, in keeping with the 
notion of forward and backward fluxes. This 
would be a necessary requirement if the two 
split-flux Jacobians would commute. It is 
hoped that by relaxing this requirement a 
stable and monotone FVS scheme could still 
be formed. We shall base our stability anal- 
ysis on the first-order (cartesian) upwind- 
differencing operator 

A.rRes = -{(F+-F+.,) 
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+ (F- + ,-F-)} (26) 

with Fourier transform 


AxA(/?, M) = - {B + (l - e~ i0 ) 

+B" (e ,/J — l) } • (27) 

For stability of the differential equation 


d U 
dt 


= AU 


(28) 


all eigenvalues of A must have a non-positive 
real part; this must be enforced for all j3 and 
A L 

In addition to stability, the split fluxes used 
in the discrete residual, (26), must yield 
monotone solutions of the state vector, U. 
Consider the steady form of (26), 


0 = A_F+ + A+F- (29) 


This equation is then linearized resulting in 
0 = B + A_U + B“A_j_U (30) 

or, for non-singular B" 

0 = BA_U + A+U, (31) 

with 

B=(B-) - 1 B + (32) 

For a monotone solution, A+U and A_U 
should have the same sign, so that the eigen- 
values of B must be all negative. As the 
eigenvalues of A and B are difficult to ob- 
tain and study analytically, we shall resort 
to their numerical evaluation for the sake of 
carrying out the stability and monotonicity 
analyses. 

First, we look at the stability properties of 
the LDFVS flux with a simple energy split- 
ting. Figure 15 summarizes the stability 
properties of the spatial operator (26) incor- 
porating the LDFVS splitting with 7 = 1.4 
and u> = 0. For 11 values of //, ranging from 



Figure 15: Maximum real part encountered 
in any of the eigenvalues of A(/5,M), versus 
//, for a range of values of //, using P = 2. 
The maximum is taken over all j3 and M. 


0 to 1, the largest positive real part encoun- 
tered among the three eigenvalues of A (/?, M) 
is plotted against v . The eigenvalues of the 
matrix A were found numerically given the 
Mach and wave numbers using a packaged 
root solver. 

For small values of fi there appears to be an 
appreciable range of v - values for which no un- 
stable eigenvalues occur. The stable range of 
v narrows down considerabty as / 1 increases, 
until only a single stable value of v remains, 
i.e. v = 0.52 for fi = 0.96. This shows that 
the diffusive fluxes at M = 0 can be reduced 
by a factor 25 without losing the possibility of 
being implemented in a stable time-marching 
scheme. Similar results, but slightly more re- 
strictive, were obtained for P = 6, so this 
avenue was not explored further. We also 
checked stability for 7=1, and found no 
qualitative difference. 

Although this result is encouraging, it is 
desirable to make the scheme strictly stable 
for standard explicit procedures. For this we 
turn to the energy splitting and perform the 
stability analysis now for fi = 1. Figure 16 
shows the maximum of all eigenvalues (over 
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Figure 16: Contours of maximum real part 
of all eigenvalues for f.i = 1. The maximum 
is taken over all (3 and M. The lowest level 
displayed is 0.002, 



Figure 17: Conical Navier-Stokes solution us- 
ing the stable splitting parameters: // = 1.0, 
v — 3/4 and lo — 3/2. 


all Mach and wave numbers) for selected val- 
ues of both v and to. One can see that the 
scheme approaches stability in a very nar- 
row trough with a minimum at v = 3/4 and 
u> = 3/2, although even there the maximum 
eigenvalue is still slightly positive. This split- 
ting is found to be stable in practice for the 
conical flow problem: for fi = 1, v = 3/4 
and io — 3/2 a converged solution can be 
obtained using a single-stage explicit (for- 
ward Euler) scheme with a Courant number 
of CFL = 0.1 after a large number (10,000) 
of iterations. Unfortunately, the stable LD- 
FVS scheme yields the worst non-monotone 
solution for the conical flow, as is shown in 
Figure 17. 

A detailed look at the eigenvalues is nec- 
essary to see if there is any hope at all of 
finding a scheme that will be both stable and 
monotone. Figure 18 shows the eigenvalues 
of B + for the original FVS with Hand’s en- 
ergy splitting and the LDFVS for \i — 1 and 
v — Lo = 0. This choice of parameters shows 
the efFect of the low diffusion mass splitting 
upon the split flux eigenvalues. This Figure 



Figure 18: Eigenvalues of the forward-flux Ja- 
cobian, B + , for the original FVS with Hand’s 
energy flux, and the LDFVS with fi = 1, 

v — u) — 0 . 
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Figure 19: Eigenvalues of the forward-flux 
Jacobian, B + , for the LDFVS with \l = 1, 

v — 3/4,o; — 3/2 and the HOPE FVS. Figure 20: Eigenvalues of the B matrix for 

the original FVS scheme using Hanel’s energy 
flux. 


shows that the eigenvalue corresponding to 
the characteristic speed u has the incorrect 
sign for M < 0, while the eigenvalue corre- 
sponding to u — a, which for M < 1 should 
be zero or slightly positive, (in the context of 
a forward flux) has the wrong sign over the 
whole Mach number range. The stable split- 
ting (// = 1 = 3/4 and o? = 3/2) and the 

HOPE splitting are compared in the next Fig- 
ure. It is seen from this Figure that the eigen- 
value corresponding to u — a still has the in- 
correct sign over a large Mach number range 
for both splittings. With this type of behav- 
ior, it is highly unlikely that a monotone solu- 
tion could generally be achieved with either 
of the non-diffusive fluxes, regardless of the 
pressure or energy splittings. This is also sup- 
ported by the results of the full monotonicity 
analysis. 

Figures 20, 21, 22, and 23 show the eigen- 
values of B over the range of Mach numbers 
|M| < 1. (These were computed by numeri- 
cally finding the inverse of B - for M < 0, or 
B + for M > 0, multiplying the inverse with 
the remaining Jacobian and finding the eigen- 
values using a root solver.) The first figure 



Figure 21: Eigenvalues of the B matrix for 
the LDFVS splitting using fi = 1, v = u = 0. 
The curves for all three eigenvalues fall nearly 
on top of each other. 
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Figure 22: Eigenvalues of the B matrix for 
the stable LDFVS splitting (ft = 1, v = 3/4 
and to — 3/2). 




- 1.00 - 0.33 0.33 1.00 

Mach 


Figure 23: Eigenvalues of the B matrix for 
the HOPE splitting. 


shows the eigenvalues for the original FVS 
scheme. As can be seen here, all of the B 
eigenvalues have the proper negative sign, as 
they should, which provides some level of con- 
fidence in the monotonicity argument. Figure 
21 shows the eigenvalues of B for the LDFVS 
with ft — 1 and v = u = 0, while Figure 22 
shows the eigenvalues for the stable LDFVS. 
Both of these schemes have all eigenvalues of 
B of the wrong sign, which indicates that a 
monotone solution, whether stable or not, is 
generally impossible and a lucky coincidence 
at best. Figure 23 shows the eigenvalues for 
HOPE. Although it does have eigenvalues of 
the wrong sign over a wide range of Mach 
numbers, there is a region near \M\ = 1 
where the correct signs are observed. This 
behavior can be traced back to Figure 19, 
which shows the split flux eigenvalues. For a 
small region near |A/| = 1 the HOPE scheme 
has the proper sign for all of the split-flux 
eigenvalues. But, as in the LDFVS schemes, 
the HOPE scheme has one eigenvalue that is 
hopelessly of the wrong sign near M — 0. 

This type of behavior cannot be avoided, 
no matter what values arc used for the pa- 
rameters in LDFVS or HOPE, including the 
exponents P and S\ The bottom line is that 
splittings that are non-difFusive for M = 0, 
by necessity must connect Van Leer’s split 
fluxes for |A/| — > 1 to central difference fluxes 
for \M\ — » 0 (see Equations (12-15)), and the 
latter are known to yield oscillatory solutions. 

7 Conclusions 

In this paper we have compared various low- 
difFusion flux formulas for use in Navier- 
Stokes computations. To compare these var- 
ious formulas, the viscous, hypersonic flow 
over a 10° cone at M = 7.95 was computed 
by solving the conical Navier- Stokes equa- 
tions in one dimension. The Harten-Lax/Roe 
scheme, incorporating a ’’smart” scalar diffu- 
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sion coefficient, presented but not tested in 

[1], was tested here and shown to give good 
results, although there is some question as 
to its robustness. Next, recent developments 
in flux- vector and hybrid flux-vector/flux- 
difference splitting (FVS/FDS) were exam- 
ined. These fluxes were shown to yield a 
correct temperature profile in the boundary 
layer, accompanied by a pressure irregular- 
ity near the boundary layer edge. In an at- 
tempt to overcome this deficiency, it was in- 
vestigated whether a pure FVS could be con- 
structed that would cause minimal dissipa- 
tion across a contact discontinuity, yet still 
be stable and monotone. We have shown 
that it is indeed possible to construct a pure 
FVS scheme with minimal mass diffusion, al- 
though this is achieved at marginal stability 
and a loss of monotonicity, even for first-order 
spatial differencing. This makes pure FVS a 
dead-end street; the latest developments in 
hybrid FVS/FDS formulas, though, indicate 
that some elements of FVS may survive and 
be used advantageously in flux functions for 
the Navier-Stokes equations. 
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